# Replication code for Divided Armies Replication and AME Plots for CoW War Types 
# projectmars V1.1 dataset 
library(readxl)
library(ggplot2) 
library(vcov)
library(sandwich)
library(margins)
library(estimatr)
projectmars <- read_excel("Desktop/projectmars.xls", na = "NA")
View(projectmars)


## Replication models with Gibler and Miller (2021) preferred specification (hereafter, "GM") 

## MIC and logged LER, GM preferred model, by CoW War Type

log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, cow_inter == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum1 <- summary(margins(log_lermean, variables = "mic_mean"))
log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean    + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, cow_intra == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum2 <- summary(margins(log_lermean, variables = "mic_mean"))
log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean     + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, cow_extra == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum3 <- summary(margins(log_lermean, variables = "mic_mean"))
log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, cow_nonstate == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum4 <- summary(margins(log_lermean, variables = "mic_mean"))
log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init +  allypowermean     + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, notincow == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum5 <- summary(margins(log_lermean, variables = "mic_mean"))

## Combine into plotting df
mic_mean.margplot <- as.data.frame(rbind(log_lermean_sum1,
                                         log_lermean_sum2,
                                         log_lermean_sum3,
                                         log_lermean_sum4,
                                         log_lermean_sum5))
mic_mean.margplot$outcome <- c("Inter-State War", #outcome labels
                               "Intra-State War",
                               "Extra-State War",
                               "Non-State War",
                               "Not in CoW")
mic_mean.margplot$outcomeorder <- seq(from = nrow(mic_mean.margplot), to = 1, by = -1)


## LER Plot by CoW War Type
ggplot(mic_mean.margplot, aes(x = reorder(outcome, outcomeorder),
                              y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Loss-Exchange Ratios") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))

## AME for all 5 battlefield performance measures, by CoW War Type 

## Inter-State War Models 

belowparity <- glm(belowparity ~ mic_mean + pol2 + pol2init + oppdemo7 + init +  allypowermean +  logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
belowparity_sum <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_mean"))
desert <- glm(desert ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
desert_sum <- summary(margins(desert, vcov = vcov_out, variables = "mic_mean"))
defect <- glm(defect ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
defect_sum <- summary(margins(defect, vcov = vcov_out, variables = "mic_mean"))
blocking <- glm(blocking ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
blocking_sum <- summary(margins(blocking, vcov = vcov_out, variables = "mic_mean"))
bpi <- lm_robust(bpi ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum <- summary(margins(bpi, variables = "mic_mean"))

## Combine into plotting df
mic_mean.margplot <- as.data.frame(rbind(belowparity_sum,
                                         desert_sum,
                                         defect_sum,
                                         blocking_sum,
                                         bpi_sum))
mic_mean.margplot$outcome <- c("LER Below Parity", #outcome labels
                               "Mass Desertion",
                               "Mass Defection",
                               "Blocking Units",
                               "BP Index")
mic_mean.margplot$outcomeorder <- seq(from = nrow(mic_mean.margplot), to = 1, by = -1)


## Inter-State War Plot
ggplot(mic_mean.margplot, aes(x = reorder(outcome, outcomeorder),
                              y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Inter-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))


## Intra-State War Models 

belowparity <- glm(belowparity ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
belowparity_sum2 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_mean"))
desert <- glm(desert ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
desert_sum2 <- summary(margins(desert, vcov = vcov_out, variables = "mic_mean"))
defect <- glm(defect ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
defect_sum2 <- summary(margins(defect, vcov = vcov_out, variables = "mic_mean"))
blocking <- glm(blocking ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
blocking_sum2 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_mean"))
bpi <- lm_robust(bpi ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum2 <- summary(margins(bpi, variables = "mic_mean"))

## Combine into plotting df
mic_mean.margplot <- as.data.frame(rbind(belowparity_sum2,
                                         desert_sum2,
                                         defect_sum2,
                                         blocking_sum2,
                                         bpi_sum2))
mic_mean.margplot$outcome <- c("LER Below Parity", #outcome labels
                               "Mass Desertion",
                               "Mass Defection",
                               "Blocking Units",
                               "BP Index")
mic_mean.margplot$outcomeorder <- seq(from = nrow(mic_mean.margplot), to = 1, by = -1)


## Intra-State War Plot
ggplot(mic_mean.margplot, aes(x = reorder(outcome, outcomeorder),
                              y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Intra-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))



## Extra-State War Models 

belowparity <- glm(belowparity ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
belowparity_sum3 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_mean"))
desert <- glm(desert ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
desert_sum3 <- summary(margins(desert, vcov = vcov_out, variables = "mic_mean"))
defect <- glm(defect ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
defect_sum3 <- summary(margins(defect, vcov = vcov_out, variables = "mic_mean"))
blocking <- glm(blocking ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
blocking_sum3 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_mean"))
bpi <- lm_robust(bpi ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum3 <- summary(margins(bpi, variables = "mic_mean"))

## Combine into plotting df
mic_mean.margplot <- as.data.frame(rbind(belowparity_sum3,
                                         desert_sum3,
                                         defect_sum3,
                                         blocking_sum3,
                                         bpi_sum3))
mic_mean.margplot$outcome <- c("LER Below Parity", #outcome labels
                               "Mass Desertion",
                               "Mass Defection",
                               "Blocking Units",
                               "BP Index")
mic_mean.margplot$outcomeorder <- seq(from = nrow(mic_mean.margplot), to = 1, by = -1)


## Extra-State War Plot
ggplot(mic_mean.margplot, aes(x = reorder(outcome, outcomeorder),
                              y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Extra-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .75, size=12.5, colour="black"))





## Non-State War Models 

belowparity <- glm(belowparity ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
belowparity_sum4 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_mean"))
desert <- glm(desert ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
desert_sum4 <- summary(margins(desert, vcov = vcov_out, variables = "mic_mean"))
defect <- glm(defect ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
defect_sum4 <- summary(margins(defect, vcov = vcov_out, variables = "mic_mean"))
blocking <- glm(blocking ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
blocking_sum4 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_mean"))
bpi <- lm_robust(bpi ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum4 <- summary(margins(bpi, variables = "mic_mean"))

## Combine into plotting df
mic_mean.margplot <- as.data.frame(rbind(belowparity_sum4,
                                         desert_sum4,
                                         defect_sum4,
                                         blocking_sum4,
                                         bpi_sum4))
mic_mean.margplot$outcome <- c("LER Below Parity", #outcome labels
                               "Mass Desertion",
                               "Mass Defection",
                               "Blocking Units",
                               "BP Index")
mic_mean.margplot$outcomeorder <- seq(from = nrow(mic_mean.margplot), to = 1, by = -1)


## Non-State War Plot
ggplot(mic_mean.margplot, aes(x = reorder(outcome, outcomeorder),
                              y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Non-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))




## Non-COW War Models 

belowparity <- glm(belowparity ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$notincow==1,]$ccode)
belowparity_sum5 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_mean"))
desert <- glm(desert ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$notincow==1,]$ccode)
desert_sum5 <- summary(margins(desert, vcov = vcov_out, variables = "mic_mean"))
defect <- glm(defect ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$notincow==1,]$ccode)
defect_sum5<- summary(margins(defect, vcov = vcov_out, variables = "mic_mean"))
blocking <- glm(blocking ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean    + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$notincow==1,]$ccode)
blocking_sum5 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_mean"))
bpi <- lm_robust(bpi ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum5 <- summary(margins(bpi, variables = "mic_mean"))

## Combine into plotting df
mic_mean.margplot <- as.data.frame(rbind(belowparity_sum5,
                                         desert_sum5,
                                         defect_sum5,
                                         blocking_sum5,
                                         bpi_sum5))
mic_mean.margplot$outcome <- c("LER Below Parity", #outcome labels
                               "Mass Desertion",
                               "Mass Defection",
                               "Blocking Units",
                               "BP Index")
mic_mean.margplot$outcomeorder <- seq(from = nrow(mic_mean.margplot), to = 1, by = -1)


## Not in COW War Plot
ggplot(mic_mean.margplot, aes(x = reorder(outcome, outcomeorder),
                              y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Not in CoW") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))



## MIC and logged LER, GM preferred model, by CoW War Type

log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum1 <- summary(margins(bpi, variables = "mic_mean"))

log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum2 <- summary(margins(bpi, variables = "mic_mean"))

log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean    + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum3 <- summary(margins(bpi, variables = "mic_mean"))

log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum4 <- summary(margins(bpi, variables = "mic_mean"))

log_lermean <- lm_robust(log_lermean ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean    + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), se_type = "stata", clusters = ccode)
summary(log_lermean)
log_lermean_sum5 <- summary(margins(bpi, variables = "mic_mean"))


## Robustness check: MIC Bands and Battlefield Performance, by CoW War Type

## Inter-State War Models 

belowparity <- glm(belowparity ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
belowparity_sum6 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_bands"))
desert <- glm(desert ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
desert_sum6 <- summary(margins(desert, vcov = vcov_out, variables = "mic_bands"))
defect <- glm(defect ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
defect_sum6 <- summary(margins(defect, vcov = vcov_out, variables = "mic_bands"))
blocking <- glm(blocking ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_inter==1,]$ccode)
blocking_sum6 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_bands"))
bpi <- lm_robust(bpi ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_inter == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum6 <- summary(margins(bpi, variables = "mic_bands"))

## Combine into plotting df
mic_bands.margplot <- as.data.frame(rbind(belowparity_sum6,
                                          desert_sum6,
                                          defect_sum6,
                                          blocking_sum6,
                                          bpi_sum6))
mic_bands.margplot$outcome <- c("LER Below Parity", #outcome labels
                                "Mass Desertion",
                                "Mass Defection",
                                "Blocking Units",
                                "BP Index")
mic_bands.margplot$outcomeorder <- seq(from = nrow(mic_bands.margplot), to = 1, by = -1)


## Inter-State War Plot
ggplot(mic_bands.margplot, aes(x = reorder(outcome, outcomeorder),
                               y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Inter-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 6)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))


## Intra-State War Models 

belowparity <- glm(belowparity ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
belowparity_sum7 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_bands"))
desert <- glm(desert ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
desert_sum7 <- summary(margins(desert, vcov = vcov_out, variables = "mic_bands"))
defect <- glm(defect ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
defect_sum7 <- summary(margins(defect, vcov = vcov_out, variables = "mic_bands"))
blocking <- glm(blocking ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_intra==1,]$ccode)
blocking_sum7 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_bands"))
bpi <- lm_robust(bpi ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_intra == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum7 <- summary(margins(bpi, variables = "mic_bands"))

## Combine into plotting df
mic_bands.margplot <- as.data.frame(rbind(belowparity_sum7,
                                          desert_sum7,
                                          defect_sum7,
                                          blocking_sum7,
                                          bpi_sum7))
mic_bands.margplot$outcome <- c("LER Below Parity", #outcome labels
                                "Mass Desertion",
                                "Mass Defection",
                                "Blocking Units",
                                "BP Index")
mic_bands.margplot$outcomeorder <- seq(from = nrow(mic_bands.margplot), to = 1, by = -1)


## Intra-State War Plot
ggplot(mic_bands.margplot, aes(x = reorder(outcome, outcomeorder),
                               y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Intra-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 6)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))



## Extra-State War Models 

belowparity <- glm(belowparity ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
belowparity_sum8 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_bands"))
desert <- glm(desert ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
desert_sum8 <- summary(margins(desert, vcov = vcov_out, variables = "mic_bands"))
defect <- glm(defect ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
defect_sum8 <- summary(margins(defect, vcov = vcov_out, variables = "mic_bands"))
blocking <- glm(blocking ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_extra==1,]$ccode)
blocking_sum8 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_bands"))
bpi <- lm_robust(bpi ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_extra == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum8 <- summary(margins(bpi, variables = "mic_bands"))

## Combine into plotting df
mic_bands.margplot <- as.data.frame(rbind(belowparity_sum8,
                                          desert_sum8,
                                          defect_sum8,
                                          blocking_sum8,
                                          bpi_sum8))
mic_bands.margplot$outcome <- c("LER Below Parity", #outcome labels
                                "Mass Desertion",
                                "Mass Defection",
                                "Blocking Units",
                                "BP Index")
mic_bands.margplot$outcomeorder <- seq(from = nrow(mic_bands.margplot), to = 1, by = -1)


## Extra-State War Plot
ggplot(mic_bands.margplot, aes(x = reorder(outcome, outcomeorder),
                               y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Extra-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 6)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))





## Non-State War Models 

belowparity <- glm(belowparity ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean  + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
belowparity_sum9 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_bands"))
desert <- glm(desert ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
desert_sum9 <- summary(margins(desert, vcov = vcov_out, variables = "mic_bands"))
defect <- glm(defect ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
defect_sum9 <- summary(margins(defect, vcov = vcov_out, variables = "mic_bands"))
blocking <- glm(blocking ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$cow_nonstate==1,]$ccode)
blocking_sum9 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_bands"))
bpi <- lm_robust(bpi ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, cow_nonstate == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum9 <- summary(margins(bpi, variables = "mic_bands"))

## Combine into plotting df
mic_bands.margplot <- as.data.frame(rbind(belowparity_sum9,
                                          desert_sum9,
                                          defect_sum9,
                                          blocking_sum9,
                                          bpi_sum9))
mic_bands.margplot$outcome <- c("LER Below Parity", #outcome labels
                                "Mass Desertion",
                                "Mass Defection",
                                "Blocking Units",
                                "BP Index")
mic_bands.margplot$outcomeorder <- seq(from = nrow(mic_bands.margplot), to = 1, by = -1)


## Non-State War Plot
ggplot(mic_bands.margplot, aes(x = reorder(outcome, outcomeorder),
                               y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Non-State War") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))




## Non-COW War Models 

belowparity <- glm(belowparity ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), family = binomial(link = "logit"))
summary(belowparity)
vcov_out <- vcovCL(belowparity, cluster = projectmars[projectmars$notincow==1,]$ccode)
belowparity_sum10 <- summary(margins(belowparity, vcov = vcov_out, variables = "mic_bands"))
desert <- glm(desert ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), family = binomial(link = "logit"))
summary(desert)
vcov_out <- vcovCL(desert, cluster = projectmars[projectmars$notincow==1,]$ccode)
desert_sum10 <- summary(margins(desert, vcov = vcov_out, variables = "mic_bands"))
defect <- glm(defect ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), family = binomial(link = "logit"))
summary(defect)
vcov_out <- vcovCL(defect, cluster = projectmars[projectmars$notincow==1,]$ccode)
defect_sum10<- summary(margins(defect, vcov = vcov_out, variables = "mic_bands"))
blocking <- glm(blocking ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean    + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow== 1), family = binomial(link = "logit"))
summary(blocking)
vcov_out <- vcovCL(blocking, cluster = projectmars[projectmars$notincow==1,]$ccode)
blocking_sum10 <- summary(margins(blocking, vcov = vcov_out, variables = "mic_bands"))
bpi <- lm_robust(bpi ~ mic_bands + pol2 + pol2init + oppdemo7 + init  + allypowermean   + logdist + standing + fullvolunteer + composite, data=subset(projectmars, notincow == 1), se_type = "stata", clusters = ccode)
summary(bpi)
bpi_sum10 <- summary(margins(bpi, variables = "mic_bands"))

## Combine into plotting df
mic_bands.margplot <- as.data.frame(rbind(belowparity_sum10,
                                          desert_sum10,
                                          defect_sum10,
                                          blocking_sum10,
                                          bpi_sum10))
mic_bands.margplot$outcome <- c("LER Below Parity", #outcome labels
                                "Mass Desertion",
                                "Mass Defection",
                                "Blocking Units",
                                "BP Index")
mic_bands.margplot$outcomeorder <- seq(from = nrow(mic_bands.margplot), to = 1, by = -1)


## Not in COW War Plot
ggplot(mic_bands.margplot, aes(x = reorder(outcome, outcomeorder),
                               y = AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Not in CoW") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8)) +
  theme_bw() +
  theme(plot.title = element_text(hjust=0.5, size=16)) + 
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=2),
        # legend.justification = c(1, 1),
        # legend.position= c(1, 1),
        legend.position = "none",
        axis.text.y = element_text(size=12.5, colour = "black"),
        axis.text.x  = element_text(angle=0, vjust=1, hjust = .5, size=12.5, colour="black"))


## Extension: Reanalysis with Clodfelter 2008-only Wars 

belowparity <- glm(belowparity ~ mic_mean + pol2 + pol2init + oppdemo7 + init + allypowermean + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, clodfelter08 == 1), family = binomial(link = "logit"))
summary(belowparity)
desert <- glm(desert ~ mic_mean + pol2 + pol2init + oppdemo7 + init + allypowermean + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, clodfelter08 == 1), family = binomial(link = "logit"))
summary(desert)
defect <- glm(defect ~ mic_mean + pol2 + pol2init + oppdemo7 + init + allypowermean + logdist + standing + fullvolunteer + composite, data=subset (ProjectMarsV1_1, clodfelter08 == 1), family = binomial(link = "logit"))
summary(defect) 
blocking <- glm(blocking ~ mic_mean + pol2 + pol2init + oppdemo7 + init + allypowermean + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, clodfelter08 == 1), family = binomial(link = "logit"))
summary(blocking) 
bpi <- lm_robust(bpi ~ mic_mean + pol2 + pol2init + oppdemo7 + init  + allypowermean + logdist + standing + fullvolunteer + composite, data=subset(ProjectMarsV1_1, clodfelter08 == 1), se_type = "stata", clusters = ccode)
summary(bpi)

